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Abstract 

We present results from lattice simulations of a monolayer graphene model at non-zero 
temperature. At low temperatures for sufficiently strong coupling the model develops 
an excitonic condensate of particle-hole pairs corresponding to an insulating phase. The 
Berezinskii-Kosterlitz-Thouless phase transition temperature is associated with the value 
of the coupling where the critical exponent 6 governing the response of the order pa- 
rameter at criticality to an external source has a value close to 15. The critical coupling 
on a lattice with temporal extent Nt = 32 (T = l/{Ntat) where at is the temporal 
lattice spacing) and spatial extent Ng = 64 is very close to infinite coupling. The value 
of the transition temperature normalized with the zero temperature fermion mass gap 
Ao is given by ^-^^ = 0.055(2). This value provides an upper bound on the transi- 
tion temperature, because simulations closer to the continuum limit where the full U{4:) 
symmetry is restored may result in an even lower value. In addition, we measured the 
helicity modulus T and the fermion thermal mass Ar(T), the latter providing evidence 
for a pseudogap phase with At > extending to arbitrarily high T. Analysis of the 
dispersion relation suggests that the Fermi velocity is not sensitive to thermal effects. 
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1 Introduction 



The impact of electron-electron interactions on the physics of graphene is an important 
focus of current study (for recent reviews see p]). There are simple arguments why an 
"independent quasiparticle" picture may not be adequate for certain properties. Firstly, 
since the carrier density of states vanishes in undoped graphene (the zero energy con- 
dition is only satisfied at two isolated "Dirac points" in the first Brillouin zone), the 
effects of screening are much less in graphene than in a conventional conductor, the 
main contribution coming from electron-hole pairs which increase the effective dielectric 
constant of the medium in a fashion entirely analogous to vacuum polarisation in QED. 
This means that the interaction between charged carriers remains Coulombic, i.e. long- 
ranged oc r~^. Secondly, the relative importance of quantum corrections, parametrised 
by the fine structure constant a, is much greater than in conventional QED, because 

2 

tteff = 4^;^/— where vp ^ sUo Fermi velocity and e the dielectric permittivity 

of the underlying substrate: hence = ~ 0(1), and its value depends on the 
substrate, taking a maximum value 2.16 for suspended graphene. 

These considerations have motivated the study of an effective (2 + l)d relativistic 
field theory with Nf fermion flavours for the low energy electronic excitations {Nf = 2 
for monolayer graphene) and an instantaneous Coulomb interaction between conserved 
charges, to be reviewed in Sec. H] below [21 El S]. For sufficiently strong coupling the 
theory describes a quantum critical point (QCP) at T = separating a semimetal phase 
in which charge carriers remain ungapped, from an insulating phase in which electron- 
hole exciton pairs condense in the ground state inducing a gap at the Dirac points. It 
is conceivable that the properties of the QCP dominate the effective description of low- 
energy charge transport in graphene irrespective of whether the semimetal or insulating 
phase is physically realised. 

Since the theory is strongly interacting, various non-perturbative approaches have 
been applied, including Monte Carlo simulation of an effective lattice field theory pos- 
tulated to belong to the same universality class at the QCP. In a series of papers, Drut 
and Lahde [S] have simulated a graphene field theory with staggered lattice fermions in 
which electrostatic degrees of freedom are formulated on a (3+l)-dimensional lattice, 
while the electron fields are restricted to a (2-|-l)-dimensional slice. Their results favour 
the scenario that suspended graphene with aes = 2.16 is an insulator. More recent sim- 
ulations with an improved fermion action support this scenario |6]. Two of us [7] have 
simulated an entirely 2+1-dimensional model which is in essence a non-covariant form 
of the Thirring model [8j, and showed that at infinite coupling for Nf < Nj^ = 4.8(2) 
graphene is an insulator, whereas for Nf > Nfc it is a semimetal. The results from simu- 
lations of the same model at finite coupling provided evidence that graphene in vacuum 
is an insulator [9] in agreement with [5l [6] . More recently, the authors of [10] presented 
preliminary results from Monte Carlo simulations of the tight-binding Hamiltonian on 
a hexagonal lattice. 

At nonzero temperature, universality arguments imply that the critical properties 
of a (d + l)-dimensional theory coincide with those of a dimensional classical spin 
model with the same symmetries. The contribution of non-zero Matsubara modes can 
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be absorbed into non-universal aspects of the transition. Consequently, fermions which 
satisfy antiperiodic boundary conditions and do not have zero modes are expected to 
decouple from the scalar sector. The validity of the dimensional reduction was confirmed 
with accuracy in Monte Carlo simulations of fermionic field theories such as the (2 + l)d 
Gross-Neveu model [H] and the (3+l)d Nambu— Jona-Lasinio (NJL) model [12] and 
strong couphng QCD [13]. 

There has been compelling experimental evidence [H] that at constant low tem- 
perature graphene undergoes a Berezinskii-Kosterlitz-Thouless (BKT) phase transition 
[15] when the intensity of an external magnetic field is varied. The authors of [T^ 
showed that the steep increase in the electrical resistance at the Dirac point as a func- 
tion of the magnetic field fitted accurately the essential scaling relation of the BKT 
scenario. The BKT transition occurs in two-dimensional systems with a U(l) symme- 
try and is driven by the unbinding of vortices, as reviewed in Sec. [31 The transition 
separates two phases, neither of which have long-range order: a low temperature spin- 
wave phase where vortices and antivortices form bound states and a high temperature 
plasma-like phase of unbound vortices and antivortices. An analytical approach based 
on solutions of self-consistent Schwinger-Dyson equations [16] predicted that the critical 
temperature is given by Tbkt = '^^(Tbkt)/'^ ~ ^o/S, where T{Tbkt) is the helicity 
modulus or stiffness of the order parameter at the transition temperature and Aq is 
the fermion mass gap at T = 0. However, care is needed since as shown in [TTj, in a 
model of graphene in which the full global symmetry is U(4) (expected for QED3 with 
Nf = 2) instead of U(l), the creation of "half- vortices" is energetically more favourable 
over the usual vortices. As a result, the critical temperature is driven to a lower value 

Tbkt = t^^{Tbkt)/^ = Tbkt/4:. 

In this paper we present results from simulations of our Thirring-like graphene 
model [H] at non-zero temperature. As we show in Sec. [2] on the lattice the remnant of 
the U(4)/U(2)®U(2) manifold in which the order parameter of the continuum theory 
assumes values in U(l); we therefore do not anticipate the existence of half vortices in 
our lattice model away from the continuum limit. 

The temperature in the simulation is given by T = 1/Ntat, where Nf is the lattice 
temporal extent and at the temporal lattice spacing. In a model with anisotropic in- 
teractions we anticipate that the temporal (a^) and spatial (a^) lattice spacings are not 
equal for arbitrary interaction coupling, i.e. the anisotropy factor a^/a^ is renormalised 
by quantum corrections governed by an action which treats time and space on a different 
footing. This has to be taken into account whenever deriving relations between physical 
quantities based on lattice observables; fortunately for the current study all quantities 
can be expressed in units of the temporal lattice spacing at. 

Furthermore, as we show in Sec. 14.11 the transition temperature in natural units is 
very low: i.e. T/Aq <^ 1. This drives the critical coupling at which the BKT phase 
transition occurs to a very strong value (close to the strong coupling limit) even when 
the temporal lattice size Nf = 32. This value of A^^^ is much larger than the values 
A^t = 6, . . . , 10 usually used in simulations of nonzero temperature QCD, and makes the 
study of the BKT scenario in graphene a computationally very difficult problem. On the 
basis of large- A^j arguments [7], we believe that at very strong couplings our Thirring-like 
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model should become similar to the instantaneous Coulomb interaction model [U [5] . 

The main goals of this work are: (i) to measure Tbkt/^o] (h) to obtain a first mea- 
surement of the helicity modulus T(T) for T > Tbkt and to compare with theoretical 
expectations; (iii) to measure the fermion mass gap Ay for T > Tbkt and to demon- 
strate that it remains nonzero even in the absence of long-range order through exciton 
condensation - this situation, which has been discussed theoretically in the context of 
the Gross-Neveu model [T8], is qualitatively similar to the pseudogap phase observed in 
the phase diagram of cuprate superconductors below optimal doping. 

The paper is organised as follows: In Sec. [2] we present both the continuum model 
and the lattice formulation used here, along with a discussion of its global symmetries 
and breaking patterns. In Sec. [3] we briefly review the classic BKT theory of the thermal 
phase transition in planar models with U(l) global symmetry, and discuss modifications 
if the global symmetry is expanded. In Sec. |4] we present our simulation results, and in 
Sec. 5 we summarise and discuss our conclusions. 



2 Formulation of the Model 

Our starting point is a model of relativistic Dirac fermions moving in 2+1 dimensions 
and interacting via an instantaneous Coulomb interaction. In Euclidean metric the 
action is [3l H [I6] : 

^/ 1 /" 

5*1 = ^ / dXo(fx{lpaJodoi^a + Vp^al-^'lpa + ^I^V'aToV'a) + ^ / dx^^d^ x{diVf , (1) 
a=l ^ ^ •' 

where e is the electron charge, vp the Fermi velocity, V the electrostatic potential, and 
the 4x4 Dirac matrices satisfy {7^, 7j,} = 25^^, = 0, . . . , 3 (note 73 does not appear 
in ([T])). For monolayer graphene the number of fermion flavours is Nf = 2. 

For sufficiently large coupling the description in terms of massless relativistic 
excitations may be disrupted by condensation of bound fermion-hole exciton pairs in 
the ground state, signalled by an order parameter (ipip) O5 with the result that a 
gap appears in the fermion spectrum, corresponding to a transition from a conductor 
to an insulator. The spontaneously broken global symmetry is l]{2Nf) generated by 
rotations of the form ip UVip, ip t— )■ '?/'f/^7375 1^^7573, with U acting on flavour in- 
dices a = 1, . . . , A^/ and \^ a 2 x 2 matrix generated by the set {1, 73, 75, ^7375}, where 
{liiilb} = V/i. The order parameter remains invariant under independent U(A^/) 
rotations generated by both 1 and ^7375, resulting in a breaking pattern 

\]{2Nf) \]{Nf) ® U(A^/). (2) 

At zero temperature, for Nj < Nf^ the model predicts a finite sequence of quantum 
critical points (QCPs) whose properties at the critical coupling el{Nf) depend on Nf [4]. 
The ground state is then an excitonic condensate for > e^. Numerical simulations of 
the lattice model described below find Nfc ~ 5 [7]. The QCP is an ultraviolet-stable 
fixed point of the renormalisation group, implying a divergent correlation length and 
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algebraic behaviour of correlation functions which in principle may be distinct from that 
of free-field theory. If the physical value of in graphene were numerically close to 
the fixed-point value, in either subcritical or supercritical regimes, then the QCP might 
dominate the behaviour of low energy charged excitations, with profound impact on the 
description of transport. Ultimately this must be settled by experiment. 

The possible relevance of a QCP has motivated the application of lattice gauge 
theory simulation techniques to the study of graphene. In this paper, we study a model 
discretised on a 2 + 1 dimensional Euclidean cubic lattice with action which for Nf = 2 
can be written in the staggered fermion formulation in the form (with bare Fermi velocity 
VF = 1) [7, 9J: 

Siatt = \Y1 Xlv^.xil + iS^,oV^)xl+il-xlv^.x{l-iS^,oV.^_o)xl-fl + m ^ xlxl + ^ 

xfii xi X 

(3) 

Here X ^^e single component Grassmann fermion fields defined on lattice sites, m an 
artificial mass gap introduced to regularise IR fluctuations on a finite system volume, and 
V a boson field, which mimics the electric potential of ([1]) in the limit g"^ ^ oo, defined 
on the links emanating from the sites in the timelike direction. The Kawamoto-Smit 
phases rjf^x = (— 1)^°'' '"^^-i are lattice analogues of the Dirac 7-matrices. Note that \4 
couples to a charge density Jqx which is the timelike component of a conserved current 
Jfj.x = ^^^[XxXx+fi + XxXx-fi]- Since V appears in Gaussian form it may be integrated out 
to yield a model of self-interacting fermions resembling the Thirring model, with a local 
interaction term of the form q^Jq^- finite the V field couples to a light, tightly- 
bound electron-hole meson which becomes massless in the limit (7^ — )■ C)0 [7| yielding 
identical dynamics to the electric potential of the gauge theory ([T]). The simulation 
results presented in Sec. Hlwere obtained not far from this limit. 

A distinct model, with an identical (2 -|- l)d fermion sector this time interacting with 
abelian lattice gauge fields defined on a (3-|-l)-dimensional lattice, has been studied by 
Drut and Lahde |5j. Their formulation is designed to reproduce the action ([1]), which 
describes a long-ranged Coulomb interaction between charges. Two comments about 
the relation between the models are worth making: 

• The fermionic sectors share the same global symmetries. In the weakly coupled 
long- wavelength limit ([3]) describes Nf = 2 four-component Dirac fermions [TU] . 

• The continuum theories modelled coincide in the strong coupling (e^,(7^ — > 00) 
and/or large- iVj limits. 

In particular, the estimate Nfc = 4.8(2) obtained using ([3]) is expected to hold for both 
models ^ . 

Next we discuss symmetry breaking in the model ([3]). In the limit m — )■ there is a 
global "chiral" symmetry 

Xx ^ exp{iaex)Xx; Xx ^ exp{iaex)Xx (4) 
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where Ex = (—1)^0+^1+^2, i^j^g lattice analogue of 75, distinguishes odd and even sublat- 
tices. For iV species of lattice fermion corresponding to N j = 2N continuum flavours, 
excitonic condensation of the form (xx) = V~^dln Z/dm 7^ is the partition func- 
tion on the Euclidean spacetime lattice) induces a spontaneous symmetry breaking of 
the form 

V{Nf/2) ® V{Nf/2) V{Nf/2). (5) 

Only in the weak-coupling continuum limit must we necessarily expect a restoration of 
the continuum breaking pattern ([2]), implying in particular that |iV| would-be Goldstone 
modes remain massive for non-zero lattice spacing [20]. At the QCP, however, weak 
coupling cannot be assumed; moreover the effective theory need not even be Lorentz 
invariant. It remains unclear, therefore, whether the enhanced symmetry of ([1]) will be 
fully restored, and a more systematic study of the discretised action as advocated in [6] 
will ultimately be needed to resolve this issue. 

Finally, we mention an important technical issue concerning the model ([3]) which does 
not apply to the gauge-theory formulation For the action ([3]) there is no symmetry 
guaranteeing transversity of the vacuum polarisation tensor (i.e. A'Tlfj^^.^ 7^ 0, where 

is the backward difference operator), resulting in an additive renormalisation of the 
coupling g"^: 

where gfi^iNj) < 00 defines the effective location of the strong coupling limit. Unitarity 
is violated for g"^ > gf^^. In refs. [TJ [21] gf^^ was identified numerically with g'^^^, defined 
by the (m- and volume-independent) location of a peak in the order parameter (xx) 
found in the broken symmetry phase. 



3 Theoretical Expectations at Nonzero Temperature 

In the excitonic phase which forms at T = for g"^ > g^, for Nj = 2 the order parameter 
(xx) = = 00^*^ spontaneously breaks a U(l) global symmetry of the action ([3]). For 
T > long-range order is forbidden by the Coleman-Mermin- Wagner theorem |22j : 
rather, we expect at low T a phase where low energy phase fluctuations are described 
by an effective Hamiltonian 

Heff oc ^(V0)* ■ (V0) ^ ^(V^)^ (7) 

where in this context the parameter T is called the helicity modulus, and correlation 
functions decay algebraically: 

lim (0(O)0t(r)) = 02(e'^(°)e-'^('')) oc r-^ (8) 

m—^O 

with critical exponent rj = T/(27rT). As temperature rises topologically non-trivial 
excitations become important. A vortex of charge q has the form (in polar coordinates 
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r, '?/') 9 = qijj, \V9\ = q/r, and energy 

E, = 7rTqHn^, (9) 

as 

where Lg is the spatial extent of the universe and the lattice spacing. Overall charge 
neutrality is thus a requirement at low T if is to remain finite. Since a vortex can be 
located at any one of {Ls/asY (dual) lattice sites, the entropy 

5 = 21n— . (10) 

The free energy F = E—TS of a |g| = 1 vortex thus changes sign at a critical temperature 

Tbkt = |t. (11) 



This is the celebrated Berezinskii-Kosterlitz-Thouless transition [15] between a low-T 
critical phase in which vortices can only exist in tightly-bound dipole pairs, and a gapped 
phase where unbound vortices form a "topological plasma" which screens the long-range 
inter-vortex interaction. 

The relation ( ITTl) remains true in a more sophisticated renormalisation group treat- 
ment [23], except that T must be replaced by its screened value T(Tbkt) exactly at the 
transition. The critical exponent rj describing correlations for T < Tbkt thus obeys 

V<Vc=^- (12) 

A related exponent S describes the response of the order parameter to a small symmetry- 
breaking explicit mass gap m via (0) oc ms . It is related to r] via the hyperscaling relation 
S = (4 — 77) /?7, yielding 

6> 6c = 15. (13) 

This picture may need modification when applied to ([1]). Aleiner et al [17] have 
performed a similar analysis for the U (2)- valued {^ip) using a Hamiltonian with inde- 
pendent moduli for U(l)- and SU(2)-valued fluctuations of the order parameter field. 
The crucial point is that the SU(2) sigma model is asymptotically free, implying that 
Tsu{2) rapidly runs to zero as high-momentum modes are integrated out, with the result 
that the U(l) effective Hamiltonian ([7]) is adequate for describing physics at large dis- 
tances. However, the richer symmetry of the order parameter permits the existence of 
a new kind of topological excitation called a half- vortex with g = ±|, whose energy is 
still given by ([2]), and which is thus much more readily formed by thermal fluctuations. 
The BKT transition temperature is accordingly modified to 

Tbkt = ^T, (14) 

with new values Vc = jq and 6c = 63. 
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4 Numerical Results 



In this section we present results from our numerical investigation of the model discussed 
in the previous section at nonzero temperature. More specifically we estimate the phys- 
ical critical temperature, detect fermion mass generation in the high temperature phase 
and study the behaviour of T at high T. In Euclidean field theory the temperature T 
is related to the time-extent Lt of the universe via T = L~[^ = {Ntat)~^ where in the 
second step a timelike lattice spacing at is specified. In general numerical simulations are 
performed with A^^ fixed, so that T is varied through variation of at{g'^). Since — ?■ 
at the QCP located at the bulk critical point g1, we deduce that in the semimetal phase 
the range < T < oo maps to the range < < g^, whereas in the insulating phase 
the same temperature range is mapped to oo > g^ > g^. In this paper we are con- 
cerned with the latter case; bearing in mind the usual convention of presenting results 
in terms on inverse coupling, and also the additive coupling renormalisation described 
in the previous section, we will therefore be working in the range g^^^^ < g^^ < g^"^. 



4.1 BKT Transition 

The first set of simulations were performed with a lattice temporal extent Nt = 16 
and spatial extents A^^ = 32,48. For these lattice volumes S'peak ~ 0.375; recall that 
the value g^^ corresponding to the infinite coupling limit has previously been identified 
with g~^^]^- However, this value of g~^^y^ is higher than the value (^pg^k ~ 0.30(2) found at 
T = ^ . Although the existence of g'^^^y- defining the effective strong coupling limit is a 
ultraviolet (UV) artifact and therefore should not depend on Nt, when Nt is comparable 
to the lattice spacing at, i.e. the UV scale becomes comparable to the IR scale, then it 
becomes difficult to disentangle the bulk and thermal transitions. In Fig. [1] we present 
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Figure 1: (color online) Exciton condensate (xx) versus m from simulations at g^^ — 0.375 on 16 x 32^ 

and 16 x 48^ lattices. 
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results for the exciton condensate (xx) versus m for Ns = 32,48 and g~'^ = 0.375. It 
appears that finite volume effects are negligible down to m = 0.00125. We then fitted 
the data at g~'^ = 0.375, 0.400 from simulations on a 16 x 32^ lattice to the scaling 
relation: 

ixx) = Cm'l'. (15) 

At the critical temperature Tbkt we expect 5 = 15. The results for the exponent (5 
and the fit qualities (x^/dof) are presented in Table [H The data and the fitted curves 
are shown in Fig. [2l The very low fit qualities and the values of b = 5.5(1), 5.1(1) for 
= 0.375 and 0.400, respectively, imply that even at fifpeak the temperature is higher 
than Tbkt: we can never go down to Tbkt in simulations with Nt = IQ. 

Table 1: Results from fits of (xx) vs m from simulations on 16 x 32^ lattices. 
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Figure 2: (color online) (xx) versus m from a 16 x 32^ lattice. 



Table 2: Results from fits of (xx) vs m from simulations on 32 x 64^ lattices. 
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Figure 3: (color online) (xx) versus m from a 32 x 64^ lattice. 



These preliminary simulations teach us that it will require very large lattices to 
identify a BKT transition. In order to approach Tbkt we tried iV^ = 32 and A^^ = 64. The 
simulations on such a large lattice at strong couplings required enormous computational 
time because the number of iterations of the conjugate gradient algorithm required for 
the inversion of the Dirac matrix kernel of ([3]) increased dramatically. For this reason it 
has not proved possible to identify a transition via singular behaviour of the susceptibility 
dixx) I oi' the specific heat as was done, say, for fermion pairing leading to long-ranged 
phase coherence in the (2 + V)d Gross-Neveu model p4|, with Tbkt/ ~ 0.5, using 
ATj = 4, AT, = 30, . . . , 150. 

Our strategy for locating Tbkt is therefore based entirely on the critical scaling 
relation f|T5|) . The data for (xx) versus m were fitted to f|T5|) for the ranges m = 
0.0025,. ..,0.010 for g'^ = 0.325, m = 0.0025, 0.0175 for g-^ = 0.350 and m = 
0.0025, 0.015 for g^"^ = 0.375. The results are presented in Table [2] and Fig. [3] shows 
the data and the fitted curves. The value oi 6 = 15.0(3) found at g~^ = 0.350 implies 
that the BKT transition occurs at this coupling. It increases to 19.1(8) at g^^ = 0.325 
which corresponds to a larger lattice spacing at and hence lower T, consistent with 
the BKT scenario. Note also that at the lowest temperature {g'"^ = 0.325) the scaling 
region shrinks as compared to higher T {g'"^ = 0.350), because as m increases the system 
crosses over to the T = scaling. The slightly increased x^/dof for g^"^ = 0.375 provides 
evidence that for g~^ > 0.350 the critical scaling based on ( IT^ is not valid because this 
coupling lies in the high temperature phase. 

In order to eliminate the lattice spacing and estimate the physical critical temperature 
at the BKT transition we measured the T = fermion mass at g'"^ = 0.350. Using point 
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Figure 4: (color online) Fermion correlator for g^^ = 0.35, m = 0.01 on a 48 x 24^ lattice. 
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Figure 5: (color online) Fermion mass gap Mf{m) versus m from simulations with g ^ ~ 0.35, 0.375 on 

a 48 X 24^ lattice. 



sources we calculated the zero momentum fermion timeslice correlator 



(16) 



where "even" refers to sites with spatial coordinate x obeying (— 1)"^^ = (—1)^^ = 1- 
This restriction improves the signal to noise ratio, and originates in the observation 
that the action ([3]) is invariant only under translations by an even number of lattice 
spacings. The simulations were performed on cold lattices with Nt = 48 and A^^ = 24 for 
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m = 0.01,0.02,0.03. In Fig. H] we present the data for Cf{t) for m = 0.01. The fermion 
correlator data were fitted to: 



Cfit) = A[exp{-Mft) - i-iyexp{-Mf{Nt-t))]. (17) 

This form assumes that the spectral density p(s) is saturated by a pole at s = Mj in both 
particle and hole branches, appropriate for zero doping. In practice this assumption is 
justified by the quality of the fit, evident in Fig. |H The minus sign between the forward 
and backward terms is due to our choice of antiperiodic boundary conditions in the 
timelike direction. The values Mf{m) extracted from fits to f|T7|) were fitted to a linear 
scaling relation Mf{m) = Aq + aim, where Aq is the mass gap. The data and the fitted 
line is shown in Fig. [51 The extrapolation to m = at g'"^ = 0.35 yields A^at = 0.57(2). 
The physical estimate for the BKT temperature is then given by: 

This result is slightly below half of the analytical prediction Tbkt/'^o ~ 1/8 obtained 
by self-consistent solution of Schwinger- Dyson equations in [16]. It is only possible to 
convert it into physical units indirectly, using the estimate Aq ~ 35 meV obtained in 
[25] by modelling the T-dependence of electrical conductivity measured in suspended 
graphene samples [26]. This yields Tbkt ~ 20 Kelvin. It should be stressed that this 
result has still to be extrapolated to the continuum limit A^t — ?■ oo, — )■ 0. Another 
factor to bear in mind once lattice discretisation artifacts disappear is that the U(4) 
global symmetry of the continuum model ([1]) will be recovered. In that de- 
scribed in Sec. [3] the critical temperature Tbkt will be smaller than the value f|T8|) by 
a factor of four, because half-vortices will become energetically favoured and dominate 
the disruption of long-range phase coherence [T7] . 

4.2 Helicity Modulus 

Next we present numerical estimates of T(T): we briefly review the method, adapted 
from [27]. The mass term in is replaced by a spatially- varying source of the form 
j exp{i9{x)ex), where the single- valued phase is defined by 

6{xi, X2) = —{nixi + n2X2). (19) 

The helicity modulus parametrises the response of the axial current J^^ = ^^^[Xx{£X)x+fi+ 
Xx{£X)x-fi], which is conserved in the limit j — )• 0: 

>(j) = T(j)W = ^(ni,n2). (20) 

To make contact with the theoretical considerations discussed above requires the extrap- 
olation j — )■ 0. Note that because V ■ J" has the same form as the kinetic energy term 
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in the action ([3]), the dimensionless variables appearing in (120|) are J'^agat, and Ta^, 
meaning that T naturally scales like a mass gap. In practice to minimise discretisation 
artifacts we choose rii = 1, n2 = 0. For technical reasons associated with the choice 
Nf = 2, the results for T presented in this paper were calculated in the "partially- 
quenched" approximation, in which equilibrated field configurations were generated us- 
ing a spatially-constant mass m, the spatially-varying source only being introduced for 
the measurement of J". 
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j 

Figure 6: (color online) T versus j from simulations with = 0.45, m = 0.00125 on 16 x 32^ and 

16 X 48^ lattices. 

Given that T is noisier than (xx) '^6 restricted our simulations to a lattice with 
Nt = 16 and were therefore only able to study high temperatures. In Fig. [6] we present 
T(j) for m = 0.00125 and g^"^ = 0.45 extracted from simulations with Lg = 32 and 
Lg = 48. It is inferred that effects due to finite are small, in contrast to results from 
the Gross- Neveu model at non-zero baryon density with T < Tbkt In order to 

extract the m = value of T for each value of j we performed linear extrapolations 
using T(m, j) = T(m = 0, j) + The results for T(m = 0, j) versus j for different 
< g'"^ corresponding to T > Tbkt are shown in Fig. [71 

Unfortunately, we don't have a model permitting a reliable extrapolation of these 
data to j — )■ 0. The data show a marked downward curvature as j — )■ and it is 
therefore plausible, bearing in mind the insensitivity to L^, that T vanishes in this 
limit, as expected for T > Tbkt (however the figure, including the point where curves 
corresponding to differing temperatures intersect at j ^ 0.125, is qualitatively very 
similar to data taken with finite Lg and fixed T < Tbkt but varying baryon density 
in the 2+ld Gross-Neveu model [2Z])- For j < 0.1 there is a clear T-dependence. For 
reference Eqn. (ITTl) predicts T(TBKT)0't = 0.040, of the same order of magnitude as T(j) 
around the "knee" seen in the data of Fig. [12] at j ~ 0.1; even though a quantitative 
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Figure 7: (color online) Chirally extrapolated T versus j for different values oi g ^ extracted from 

simulations on a 16 x 32^ lattice. 



description is still lacking, therefore, the signal for T is broadly consistent with the BKT 
scenario outlined in Sec. El 

4.3 Quasiparticle Thermal Mass and Dispersion Relation 

Next, we calculated the fermion thermal mass in the high temperature region from 
simulations on 16 x 32^ lattices. Once again, the fact that the fermion correlator has a 
smaller signal-to-noise ratio than the order parameter (xx) forces us to work on smaller 
volumes. Now, at T > fermions can acquire a non-zero thermal mass even in the 
absence of spontaneous symmetry breaking. For a weakly-coupled theory, this is simply 
the Debye screening mass ~ gT, but in a strongly-coupled theory where dynamical 
mass generation at T = results from spontaneous symmetry breaking, it is better to 
draw analogies with the "pseudogap" phase thought to form in cuprate superconductors 
at strong coupling or low carrier density [l8j. Once again, we write the pairing field as 
XX = 0oe*^. For a temperature range Tbkt < T < T*, the pseudogap phase arises due 
to the "local" gap modulus (po, neutral under f/(l) rotations, remaining nonzero, while 
the phase 6 fluctuates violently, precluding both a non-zero order parameter and also the 
long-ranged phase coherence signalled by a non-vanishing helicity modulus. In Ref. [18] 
the temperature T* in the {2 + l)d Gross-Neveu model is predicted to coincide with the 
estimate Ao/2 In 2 given by mean field theory, and the difference T*—Tbkt — {Nf In 2)~^. 
The existence of the pseudogap phase at non-zero temperature was demonstrated in 
numerical simulations of Gross-Neveu models with U{1) [21] and SU{2) x SU{2) [28] 
chiral symmetries, and analytically in the Ad NJL model [29] . 

In Fig. [8] we show the fermion timeslice correlator Cj{t) for g'"^ = 0.45, 0.50, 55 and 
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13 = q'.5Q 




Figure 8: (color online) Fermion correlator with m = 0.00125 and g ^ ~ 0.45, 0.50, 0.55 on a 16 x 32^ 
lattice. The curves result from fits to data with t odd. 



m = 0.00125. We fitted the data for odd timeslices only to 

= A[exp(-Mjt) + eM-{Nt - t)Mjt)] 



(21) 



Tfie small values of Cj{t) observed on even timeslices signals a manifest chiral symmetry 
which is broken only explicitly by the fermion bare mass term. The U{1)^ symmetry (jl]) 
of staggered fermions implies that the only nonvanishing elements of the propagator are 
Cfeo and Cfoe, where the e/o subscripts denote sites with = ±1- 
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Figure 9: (color online) Fermion thermal mass Mj versus m for various g ^ extracted from simulations 

on a 16 X 32^ lattice. 
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Figure 10: (color online) Chirally extrapolated thermal mass At versus g~'^ extracted from simulations 

on a 16 X 32^ lattice. 

In Fig. Owe present the results for Mj versus m extrapolated with a linear function 
Mj{m) = At + a^m to the chiral limit. Fig. [TOl shows At versus g As g ^ increases 
the lattice spacing decreases and at the bulk critical coupling g'"^ = 0.609(2) at = Qg = 
[9j, implying T — )■ oo. It is clear from Fig. [10] that At remains of the same order of 
magnitude as Aq for a significant extent of the high temperature phase T > Tbkt, 
lending strong support to the pseudogap scenario with T* > Tbkt- 

The fermion energy as a function of momentum is accessed via analysis of the Eu- 
clidean timeslice propagator Cf{p,t) defined by 

Cjip, t)=J2 (X{0, 0)x(x, t))e-'^^-', (22) 

X even 

where the components of momentum p take values 27m/ Lg, with n = 0, 1, . . . , The 
energy E[p) is then extracted by a fit of the form 

Cf{p, t) = 5(e-^* + e-^(^*-*)), (23) 

where again only data with t odd were used. We measured E{p) for p= (pi, 0) on 16 x 32^ 
in the high temperature phase. To proceed we parametrise the dispersion relation using 



E{p) = Asmh-\^Jsm'^p + M"^), (24) 

where for A = 1 and M = m the exact result for non-interacting lattice fermions 
is recovered. Sample fits to (IMl) at m = 0.005 are shown in Fig. [TTl The dispersion 
flattens out to have zero slope at the effective Brillouin zone edge aX p = ^; this flattening 
is a discretisation artifact with no physical signiflcance. For small M we can interpret 
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£"(0) = MJ ^ AM as the quasiparticle mass (or gap), and for small p in the limit M — )■ 
then dE/dp y4 is the renormalized Fermi velocity v'^^j^at/ag at nonzero temperature, 
where we have restored explicit factors of lattice spacing. Without further information 
we are unable to distinguish between renormalization of the physical Fermi velocity and 
that of the cutoff anisotropy due to quantum corrections (this point was not realised in 
[9]), but note that the latter must be T- independent. Results for A as a function of m 
are shown in Fig. [121 Despite some noise in the data the parameter A, and hence Vp^^, is 
both m- and (^"^-independent taking a numerical value ~ 0.65, which is very close to the 
value A ^ 0.7 reported in [9] at T = 0. This implies that the principal physical effect of 
the hot medium is to generate a nonzero thermal mass, rather than to renormalize the 
Fermi velocity. A similar effect was observed in nonzero T simulations of the (2 + l)d 
Gross- Neveu model with an SU{2) (g) SU{2) chiral symmetry [28] . 
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Figure 11: (color online) Quasiparticle dispersion relation E{p) as measured on a 16 x 32^ lattice with 

TO = 0.005. 



5 Summary and Conclusion 

The main result of this first, exploratory study of thermal effects in the insulating 
phase of the graphene effective theory ([1]) with Nf = 2, via numerical simulation of its 
discrete avatar (|3]), is the determination of the critical temperature for vortex unbinding 
Tbkt/^o ~ 0.06. This value is considerably smaller than the ratio found in the Gross- 
Neveu model {Tbkt/^o ~ 0.5) |24], underlining the point that different four-point Fermi 
interactions yield distinct dynamics in (2 + l)d, and that perturbative approaches such 
as the 1/A^/ expansion are unlikely to be accurate for graphene [8]. It also implies 
that study of the BKT transition in this system is a numerically challenging problem, 
requiring large lattice volumes in order to resolve the large separation of scales. With 
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Figure 12: (color online) The fitted parameter A vs m for various values of g^^. 

the resources at our disposal we have been able to work with Nt = 32, which has enabled 
an estimate of Tbkt via the critical scaling (fT3|) of the order parameter with external 
mass source and identification of the exponent 5, but not yet, it must be stressed, via 
direct observation of singular behaviour in any thermodynamic observable. That said, 
it is noteworthy that our value f llSp is not too far removed from predictions made using 
Schwinger- Dyson equations [16j. 

Two major caveats must be noted. First, predictions made using the discrete model 
(|3]) can strictly only be applicable in the continuum limit; we therefore need to explore 
the limit \ to control the inevitable discretisation artifacts, which may scale 
with non-trivial powers of a^, at as the QCP is approached. Unfortunately in practical 
terms this requires the limit Nt ^ oo. Secondly, as noted earlier, it is argued that in 
the continuum limit the global symmetry of the effective graphene Lagrangian enlarges 
from U(1)(S)U(1) to U(4), implying the existence of half-vortex topological excitations, 
which exhibit an unbinding transition at a still lower temperature Tbkt = Tbkt/ 4: [TT] - 
Since our estimate of the critical temperature assumes the orthodox BKT scenario, we 
are unable to comment further on this possibility. Resolving this question will probably 
require a more refined lattice fermion discretisation, as advocated in p]. 

We have also presented results for the helicity modulus T as a function of the source 
strength j introduced to induce a circulating supercurrent in our system. The numerical 
challenge has so far restricted our study to the region T > Tbkt, but the magnitude of 
T(j) observed is consistent with the expectations of the conventional BKT scenario. We 
are unaware of any effective model enabling a controlled j ^ extrapolation on finite 
systems. 

Finally, the calculation of the quasiparticle propagator presented in Sec. 14.31 reveals 
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the persistence of a gap At ^ Aq for temperatures T > Tbkt, despite the fact that 
the form of the correlators shown in Fig. [8] is characteristic of propagation through a 
chirally-symmetric medium. As argued in p4j, in this phase the fermion flips chirality, 
permitting propagation at speeds v < vp, by constantly exchanging massless bosonic 
quanta with the medium: this is signalled by the spectral density function p{s) being 
modified from a simple pole on the mass shell to a branch cut above the threshold at 
s = Ay. The situation qualitatively resembles the discussion of the pseudogap phase 
in cuprates given in |T8]. In addition, the analysis of the fermion dispersion relation 
for T > Tbkt showed that the main effect of the hot medium is to generate a non-zero 
thermal quasiparticle mass rather than to renormalize the T = physical Fermi velocity. 
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